Numerical study of discrete models in the class of the nonlinear molecular beam 

epitaxy equation 
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We study numerically some discrete growth models belonging to the class of the nonlinear molecular 
beam epitaxy equation, or Villain-Lai-Das Sarma (VLDS) equation. The conserved restricted solid- 
on-solid model (CRSOS) with maximum heights differences AHmax = 1 and AH ma x = 2 was 
analyzed in substrate dimensions d = 1 and d — 2. The Das Sarma and Tamborenea (DT) model 
and a competitive model involving random deposition and CRSOS deposition were studied in d = 1. 
For the CRSOS model with AH max = 1 we obtain the more accurate estimates of scaling exponents 
in d = 1: roughness exponent a — 0.94 ± 0.02 and dynamical exponent z — 2.88 ± 0.04. These 
estimates are significantly below the values of one-loop renormalization for the VLDS theory, which 
confirms Janssen's proposal of the existence of higher order corrections. The roughness exponent in 
d = 2 is very near the one- loop result a = 2/3, in agreement with previous works. The moments 
W n of orders n = 2, 3, 4 of the heights distribution were calculated for all models and the skewness 
S ee W 3 /W2 3/2 and the kurtosis Q = W4/W2 2 - 3 were estimated. At the steady states, the CRSOS 
models and the competitive model have nearly the same values of S and Q in d = 1, which suggests 
that these amplitude ratios are universal in the VLDS class. The estimates for the DT model are 
different, possibly due to their typically long crossover to asymptotic values. Results for the CRSOS 
models in d — 2 also suggest that those quantities are universal. 

PACS numbers: 05.40.-a, 05.50.+q , 81.15.Aa 

Keywords: deposition models; thin films; molecular beam epitaxy; interface growth; universality 
classes; scaling exponents. 



I. INTRODUCTION 

Surface and interface growth processes are subjects of 
great interest for the perspective of applications to thin 
films and multilayers growth and, from the theoretical 
point of view, for their important role in Non-Equilibrium 
Statistical Mechanics [1,2]. Frequently those processes 
are described by discrete models which represent the ba- 
sic growth mechanisms by simple stochastic rules, such 
as aggregation and diffusion, and neglect details of the 
microscopic interactions. On the other hand, continuous 
theories are successful at representing those processes in 
the hydrodynamic limit. They predict the scaling expo- 
nents of many discrete models, which are consequently 
grouped in a small number of universality classes. 

Growth by molecular beam epitaxy [MBE), which is 
one of the most important techniques to produce high 
quality films with smooth surfaces, motivated the pro- 
posal of many discrete and continuous models. The dy- 
namics during MBE deposition is dominated by diffusion 
processes, which led to the proposal of an important theo- 
retical model, the Villain-Lai-Das Sarma (VLDS) growth 
equation [3,4] 



frequently called nonlinear molecular beam epitaxy equa- 
tion or conserved Kardar-Parisi-Zhang equation [1,5], 

The most important geometrical quantity to character- 
ize the surface of the deposit grown by such processes is 
the interface width. It is defined as the root mean square 
fluctuation of the average height 



{h-h)' 



1/2 



For short times, it scales as 



(2) 



(3) 



where (3 is called the growth exponent. For long times, 
in the steady state, the interface width saturates at 



£,sat ~ L a , 



(4) 



where a is called the roughness exponent. The crossover 
time from the growth regime to the steady state scales 
with L with the dynamical exponent 



a//3. 



(5) 



^ = v ^h + A 4 V 2 (V/j) 2 + v (x,t), 
at 



(1) 



where h(x, t) is the height at position x and time t in 
a c?-dimensional substrate, 1^4 and A4 are constants and 
r] is a Gaussian (nonconservative) noise. Eq. (1) is also 



For the VLDS theory, a one-loop dynamical 
renormalization-group (DRG) calculation [3,4] led to a = 
(4 - d)/3, z = (8 + G0/3 and (3 = (4 - d)/(8 + d) below 
the upper critical dimension d c = 4. See also the recent 
work of Katzav [6] , based on a self-consistent expansion 
approach, which also obtains these estimates. Some au- 
thors assumed the one-loop values to be exact in all or- 
ders, but Janssen [7] recently claimed that this conclusion 
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was derived from an ill-defined transformation and, con- 
sequently, there would be higher order corrections. From 
a two-loop calculation, he obtained small negative correc- 
tions to a and z in all dimensions [7]. Numerical studies 
of some discrete models which belong to the VLDS class 
in the continuum limit (large lattices, long times) were 
not able to solve this controversy. In d = 1, numerical 
work on a conserved restricted solid-on-solid model (to be 
defined below) systematically suggest a < 1 [8,9], but the 
error bars are large and, consequently, the authors still 
suggest the validity of the one- loop result. In d = 2 and 
higher dimensions [10], numerical results indicated that 
possible corrections to the one-loop result were smaller 
than the two- loops estimates of Janssen [7] . 

Another important question is motivated by recent re- 
sults on discrete models belonging to the Kardar-Parisi- 
Zhang (KPZ) class in d — 2. The KPZ growth equation 
includes second order linear and nonlinear terms which 
are more relevant than those in the VLDS equation (Eq. 
1) in the hydrodynamic limit [5,1]. Works on discrete 
KPZ models showed that the steady state values of the 
moments of the height distribution, 



W n = {(h-h) n ), 



obey power-counting, i. e. they scale as 



(6) 



(7) 



(note that Wi = £ 2 ). Moreover, estimates of the skew- 
ness 



S = 



W 3 



and of the kurtosis 



Q 



w 2 



W2 2 



3/2 



(8) 



(9) 



of the KPZ models indicated that the amplitude ratios of 
the moments W n (such as S and Q) are universal [11-13]. 
It seems that no previous work has considered these ques- 
tions in models belonging to the VLDS class, possibly due 
to the large times involved in their simulations (the dy- 
namical exponent is nearly the double of the KPZ value). 
Besides the theoretical relevance of those questions, ad- 
ditional motivation for their analysis is the fact that the 
amplitude ratios can be measured with much higher ac- 
curacy than the scaling exponents and may eventually 
help one to infer the universality class of an experimen- 
tal growth process. 

There is a small number of discrete models belonging 
to the VLDS class in the continuum limit. The discrete 
model proposed by Das Sarma and Tamborcnca (DT 
model) [14] is an example of a MBE-motivated model 
which falls in that class in d = 1, although there is evi- 
dence that its class in d = 2 is different [15,16]. On the 
other hand, the so-called conserved restricted-solid-on- 
solid (CRSOS) models, first proposed by Kim et al [8], is 



expected to belong to the VLDS class in all dimensions. 
This was already proved analytically in d = 1 [17-19]. In 
the CRSOS models, the difference of the heights of neigh- 
boring columns are always smaller than a certain value 
AH max , similarly to the RSOS model of Kim and Koster- 
litz [20,21]. However, in the Kim-Kosterlitz model, if the 
aggregation at the column of incidence does not satisfy 
that condition, then the aggregation attempt is rejected 
(consequently, the model is in the KPZ class). On the 
other hand, in the CRSOS model, the incident particle 
migrates to the nearest column at which the height dif- 
ference constraint is satisfied after aggregation. Thus, all 
deposition attempts are successful in the CRSOS model. 

Here, we will study numerically a modified version of 
the CRSOS model in d = 1 and d — 2, with two differ- 
ent values of AH max , the DT model in d = 1, simulated 
with noise-reduction methods, and a competitive model 
involving CRSOS and random deposition in d = 1. All 
these models belong to the VLDS class. We will perform 
systematic extrapolations of effective (roughness and dy- 
namical) exponents for the CRSOS model in d — 1 and 
d = 2. The asymptotic exponents in d = 1 are clearly dif- 
ferent from the one-loop DRG values and the sign of the 
deviations are in qualitative agreement with Janssen's re- 
sults [7]. In d = 2, possible corrections in the exponent 
a are smaller than the two-loop corrections calculated in 
that work, confirming other authors' conclusions. It will 
also be shown that the moments of the heights distribu- 
tion obey power-counting (Eq. 7) in d = 1 and d = 2, 
similarly to KPZ, and that the skewness and the kurto- 
sis for different versions of the CRSOS model (different 
A.H max ) and for the competitive model have nearly the 
same values. These estimates differ from those of the DT 
model in d = 1, but universality of amplitude ratios in 
the VLDS class cannot be discarded due to the typical 
long crossovers of the DT model. 

The rest of this paper is organized as follows. In Sec. 
II we present the stochastic rules of the CRSOS and DT 
models and give information on the simulation proce- 
dure. In Sec. Ill, we calculate the scaling exponents 
of the VLDS class in one-dimensional substrates. In Sec. 
IV, we calculate the scaling exponents in two-dimensional 
substrates. In Sec. V, we compare the asymptotic ampli- 
tude ratios of all models in d — 1 and d = 2. In Sec. VI 
we summarize our results and present our conclusions. 



II. MODELS AND SIMULATION PROCEDURE 

The rules for choosing the aggregation point in our 
version of the CRSOS model are slightly different from 
the original ones. The present version was introduced in 
Ref. [22] as a model for amorphous carbon-nitrogen films 
growth, but only small lattices were analyzed there and, 
consequently, reliable estimates of scaling exponents were 
not obtained. 

At any time, all pairs of neighboring columns are re- 
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stricted to obey the condition Ah < AH max , where Ah 
is the difference in the columns' heights and AH max is 
fixed. The deposition attempt begins with the random 
choice of one substrate column i. If the above condition is 
satisfied after aggregation of a new particle at the top of 
column i, then the aggregation takes place at that posi- 
tion. Otherwise, a nearest neighbor column is randomly 
chosen (independently of its height) and the same test 
is performed. This process is continued until a column 
is chosen in which the new particle can be permanently 
deposited. Here, the cases AH max = 1 and AH max = 2 
will be analyzed. 

In the original version of the CRSOS model [8], the 
aggregation takes place at the nearest column in which 
the condition on heights differences is satisfied, but in 
our version the incident particle performs a random walk 
along the substrate direction(s) while it searches for the 
aggregation point. The original model was proved to 
belong to the VLDS class in d = 1 by different meth- 
ods [17-19] and the coefficients of the VLDS equation 
were explicitly calculated for AH max = 1 [18,19]. Since 
our version does not change any symmetry of the orig- 
inal CRSOS model, it is also expected to be in that 
class. Notice, for instance, that there is no upward or 
downward current in our model due to the mechanism of 
random walks for choosing the aggregation position (the 
random steps do not depend on the relative heights of 
the columns). It implies that the coefficient of the sec- 
ond order height derivative of the growth equation (not 
shown in Eq. 1) is exactly zero, the VLDS equation be- 
ing the most plausible continuum description - see e. g. 
the discussion in Ref. [23]. 

We will also study the DT model in d = 1. In this 
model, the incident particle sticks at the top of the ran- 
domly chosen column i if it has one or two lateral neigh- 
bors at that position (a kink site or a valley, respectively) . 
Otherwise, the neighboring columns (at the right and the 
left sides in d = 1) are consulted. If the top position of 
only one of these columns is a kink site or a valley, then 
the incident particle aggregates at that point. If no neigh- 
boring column satisfies that condition, then the particle 
sticks at the top of column i. Finally, if both neighbor- 
ing columns satisfy that condition, then one of them is 
randomly chosen. 

In our simulations of the DT model, we used the noise 
reduction technique adopted in Ref. [24]. The noise re- 
duction factor m is the number of attempts at a site for 
an actual aggregation process to occur [25,26]. Here, the 
value to = 10 will be considered because it provided ac- 
curate estimates of scaling exponents in Ref. [24] from 
simulations in relatively small systems. On the other 
hand, the data for the original DT model present huge 
finite-size corrections (see e. g. Ref. [27]). 

In order to improve our discussion on the universality 
of amplitude ratios (Sec. V), we also simulated a compet- 
itive model in which the aggregation of the incident parti- 
cle may follow two different rules: with probability p, the 
particle aggregates at the top of the column of incidence, 



such as in the random deposition (RD) model [1]; other- 
wise (probability 1—p), it diffuses until finding a column 
i in which the condition hi — hj < AH max is satisfied for 
all nearest neighbors j after aggregation. Thus, the latter 
aggregation mechanism works for preserving the columns 
heights' constraint of the CRSOS model. Extending pre- 
vious conclusions on other competitive models [28,29], it 
is expected that this model is described asymptotically by 
the VLDS equation, similarly to the pure CRSOS model, 
but the coefficients U4 and A4 of the corresponding con- 
tinuous equation (Eq. 1) are expected to depend on p. 
In this paper, we will simulate the model with p = 0.25 
(p = is the pure CRSOS model). 

The above models were simulated in d = 1 in lattices 
of lengths ranging from L = 16 to L = 1024 for the CR- 
SOS model with AH max — 1 and AH max = 2, from 
L = 16 to L = 256 for the DT model and from L = 16 
to L = 512 for the competitive model. For the CR- 
SOS models, the number of realizations up to the steady 
state was typically 10 4 for the smallest lattices and nearly 
500 for the largest lattices. The same applies to the DT 
model, but notice that the largest length in that case 
was just L — 256. In d = 2, the CRSOS model with 
AH max — 1 was simulated in lattices of lengths ranging 
from L = 16 to L = 256, and with AH max = 2 only until 
L = 128. Whenever the number of realizations up to the 
steady state was smaller than 10 4 , a larger number of re- 
alizations covering the growth and the crossover regions 
was generated. This allowed the calculation of crossover 
times (see below) with good accuracy in d = 1. 

The calculation of the moments of the height distribu- 
tion at the steady states, W n (Eq. 6), followed the same 
lines described in Ref. [13]. In order to estimate dynam- 
ical exponents, we used a recently proposed method to 
calculate a characteristic time t$ which is proportional to 
the time of relaxation to the steady state [30] . For fixed 
L, after calculating the saturation width £ S at{L), t is 
defined through 

Z(L,t ) = kUt(L), (10) 

with a constant k < 1. From the Family- Vicsek rela- 
tion [31], it is expected that [30] 

r ~ L z . (11) 

Here, we estimated To with k ranging from k = 0.4 to 
k = 0.7. Since the exponent z is large, the characteristic 
times To increase very fast with L. Consequently, for 
large k, the accuracy of t is low in large lattices. On 
the other hand, for small k, the times To in small lattices 
are also very small (near To = 1) and, consequently, there 
are effects of the initial flat substrate. This is the reason 
why we chose a restricted range of k to analyze our data. 
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III. SCALING EXPONENTS IN 
ONE-DIMENSIONAL SUBSTRATES 

In order to estimate the roughness exponent from the 
interface width £, the first step is to calculate the effective 
exponents 



a 



(L,i) 



ln[Ut(L)/Ut (L/i)] 
Ini 



(12) 



for fixed i. It is expected that a^,i) —> a f° r any choice 
of i. 
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the effective exponents. On the other hand, estimating 
the asymptotic a is possible because there is no evidence 
of an upward curvature of those plots for large L. 
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FIG. 2. Effective roughness exponents (a) Q(l,2) and (b) 
a (L,4) versus 1/L 1/2 for the 1 + 1-dimensional CRSOS model 
with AH max = 2. Error bars are shown only when they are 
larger than the size of the data points. 



FIG. 1. Effective roughness exponents (a) ot(L,i) and (b) 
a (L,4) versus inverse lattice length for the 1 + 1-dimensional 
CRSOS model with AH max — 1. Error bars are shown only 
when they are larger than the size of the data points. 

In Figs, la and lb we show ct(L,2) and a (L,4) versus 
1/L, respectively, for the CRSOS model with AH max = 
1. The evolution of the data suggests that ct(L,i) con- 
verges to 0.91 < a < 0.94, accounting for the error bars 
and reasonable finite-size corrections. 

The type of plot in Figs, la and lb is suitable to fit 
the data to the scaling form 

a {L}i) «a + AL- A , (13) 

with A constant, if the correct variable L~ A is used in 
the abscissa (A = 1 was tested in Figs, la and lb). 
In its turn, Eq. (13) is a consequence of a scaling rela- 
tion £ sat L a (a + ciiL~ A ), with ao and a\ constants, 
which includes a sub-dominant term in addition to the 
dominant one in Eq. (4). However, no variable of the 
form L~ A provided a reasonable linear fit in the range 
of lattice size analyzed there. Thus, A = 1 was used in 
Figs, la and lb just to illustrate the L-dependence of 



The data for the CRSOS model with AH max = 2 was 
analyzed along the same lines. In Figs. 2a and 2b we 
show ot(L,2) an d ol(l,4) versus 1/L 1 / 2 , respectively. The 
variable in the abscissa of Figs. 2a and 2b was chosen to 
provide a good linear fit of the ot(L,i) data - see dotted 
line in Fig. 2b. These results suggest stronger finite-size 
corrections for a^,i) when compared to the model with 
AH max — 1. The corresponding asymptotic estimates 
are in the range 0.92 < a < 0.97, also accounting for 
the error bars. However, since these error bars are larger 
than those for AH max = 1, it is possible that the true 
asymptotic regime was not attained yet and that the true 
leading corrections are different. Anyway, those results 
still suggest that a < 1 in the L — > oo limit. 

Alternatively, we will analyze our data assuming the 
presence of a constant term as the sub-leading correction 
to the scaling of £ sa t 2 : 

Ut 2 =ei+AL 2a . (14) 

(since a ~ 1, it corresponds asymptotically to A ^ 2 in 
Eq. 13). is called intrinsic width and is frequently 
associated to large local slopes in discrete KPZ mod- 
els [25,26,13]. Effective exponents which cancel the 
contribution of £f may be defined as 



4 



L 2 



_ 1 In [fsat (2L) - g 2 at (£)] / fc 2 at (L) - (£/2)] 



In 2 



(15) 
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FIG. 3. Effective roughness exponents (accounting for 
the intrinsic width) versus 1/L for 1 + 1-dimensional CRSOS 
models with (a) AH ma x = 1 and (b) AH ma x = 2. 



In Figs. 3a and 3b we show ot'p versus 1/L for the 
CRSOS model with AH max — 1 and AH max — 2, re- 
spectively. Here, the variable 1/L in the abscissa was 
also not chosen to perform data extrapolation. The ef- 
fective exponents vary within narrow ranges (0.89 to 0.94 
for AH max = 1, 0.90 to 0.96 for AH max = 2), even in- 
cluding their error bars. Consequently, any variable in 
the form L~ A (0.5 < A < 2) leads to nearly the same 
extrapolated value of a. The data for AH max = 1 are 
more accurate and suggests 0.90 < a < 0.95, which is 
consistent with the previous analysis. The results for 
AH max — 2 confirm the trend to a < 1, although the 
uncertainties are larger. 

Assuming the power-counting property (Eq. 7) of the 
moments of the width distribution (to be discussed in 
detail in Sec. V), we may also use higher moments to 
estimate a. The effective exponents obtained from W3 
have large fluctuations, but those obtained from W4 be- 
have similarly to the ones obtained from the interface 
width. They are defined as 



a 



(4) _ ln[W 4 . sat (L)/W 4 , sat (L/i)} 
(LA ~ lni 



(16) 



where W^ sa t (L) are the fourth moments calculated at 
the steady states. 
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FIG. 4. Effective roughness exponents 
tained from the fourth moment W4) 



(L,2) ( ob " 

versus 1/L 1 / 2 for 
1 + 1-dimensional CRSOS models with (a) AH max — 1 and 
(b) AH max = 2. Error bars are shown only when they are 
larger than the size of the data points. 



In Figs. 4a and 4b we show cc^ 2 ) versus l/i 1 / 2 for 
the CRSOS models with AH max = 1 and AH max = 2, 
respectively. The variable in the abscissa of Figs. 4a and 
4b was also chosen to illustrate the behavior of the data 
for large L and not to fit the data to a certain scaling 
form. The downward curvature of the plots for large L 
also suggest a < 1 . The maximum and minimum reason- 
able limits that can be inferred from the evolution of the 
data for AH max = 1 give 0.92 < a < 0.96. The accuracy 
of the estimate for AH max = 2 is lower, as before. 

The intersection of at least two of the above estimates 
for AH max — 1, obtained from the scaling of different 
quantities and assuming different forms of finite-size cor- 
rections, provides a final estimate a — 0.94 ± 0.02. As 
will be discussed below, results for the DT model do not 
improve those obtained with the CRSOS model. 

In Figs. 5a and 5b we show the effective exponents 
Q^l.2) an d a (L2) f° r ^ nc n °isc-reduced DT model, also 
as a function of 1/L 1 / 2 . They are larger than a = 1 
and systematically increase with L. However, from all 
previous theoretical work and the above numerical data 
for the CRSOS models, there is no reason to expect a > 1 
in the VLDS class. Consequently, extrapolation of those 
data will not give reliable information for the discussion 
on the exponents of the VLDS theory in 1 + 1 dimensions. 
Instead, it is expected that the effective exponents for the 
noise-reduced DT model (Figs. 5a and 5b) will eventually 
begin to decrease with L, possibly for much larger L. 
Such decrease ola^ L 2 ) is actually observed in the original 
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DT model (without noise reduction), in the same range 
of lattice lengths analyzed here [27]. Also recall that, 
as shown in Ref. [27], the data for original DT model 
also present huge finite-size effects and cannot be used to 
obtain reliable estimates of VLDS exponents. 




1/L 1 / 2 

(a) 



1/L 1 / 2 

(b) 



FIG. 5. Effective roughness exponents (a) c*(l,2) (obtained 
from the interface width) and (b) 2 ) (obtained from Wa) 

versus 1/L 1 ^ 2 for the 1 + 1-dimensional DT model. Error bars 
are shown only when they are larger than the size of the data 
points. 



No improvement of the results in Figs. 5a and 5b is 
obtained by considering the contribution of the intrinsic 
width (Eqs. 14 and 15). 

There are other two points concerning our results for 
the DT model that deserve some comments. The first 
one is the comparison with results of Punyindu and Das 
Sarma in Ref. [24], who obtained awl with noise reduc- 
tion in lattice lengths L < 60. Our effective exponents 
for the smallest lattices (16 < L < 64) correspond to 
two data points at the left sides (larger 1/L) of Figs. 5a 
and 5b and those exponents are also near a = 1. Conse- 
quently, our estimates are consistent with those of Ref. 
[24]. On the other hand, we conclude that the noise- 
reduction scheme works properly only in a special range 
of lattice lengths, since its application to larger lattices 
(L = 128 and L = 256 in Figs. 5a and 5b) led to ef- 
fective exponents larger than 1, indicating much more 
complicated finite-size behavior. 

The other important point is related to the large er- 
ror bars, particularly for L = 256. One of the reasons is 
certainly the relatively small number of realizations for 
the largest lengths (see Sec. II). However, the surfaces 
generated by the DT model in d = 1 present grooves 
which may survive during long times. These structures 



largely increase the interface width of some realizations 
(see Ref. [32]) and, consequently, have remarkable influ- 
ence on the fluctuations of that quantity when averaged 
over various realizations. However, note that this insta- 
bility is controlled in the DT model, i. c. the depths of 
the grooves do not diverge as time increases, contrary to 
other discretized growth models which show true insta- 
bilities when pillars or grooves are formed [32,33]. 

Now we turn to the calculation of the dynamical expo- 
nent. 

Effective dynamical exponents are defined as 



In [t (£) /t (L/i)} 



(17) 



so that zl — > z as t — > oo. The error bars of to are 
larger than those of £ and the uncertainties are enlarged 
in the calculation of effective exponents for small values 
of i (Eq. 17), then we will work only with i = 4. 
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FIG. 6. Effective dynamical exponents 2(l,4) versus 1/L 



for the 1 + 1-dimensional CRSOS model with AH Tl 



= 1. 



Small horizontal shifts of the data points were used to avoid 
their superposition. Error bars (not shown) are smaller than 
Az = 0.02 (of this order for the largest L). 



In Fig. 6 we show z<l,4) versus 1/L for the CRSOS 
model with AH max = 1, with to calculated using four 
different values of k in Eq. (10) (0.4 < k < 0.7). The 
data for different k clearly converge to the same region, 
providing an asymptotic estimate z = 2.88 ± 0.04. This 
final estimate also accounts for the error bars (not shown 
in Fig. 6), which are near Az = 0.02 for the largest values 
of L. Again it is clear that the value z — 3 of one-loop 
renormalization is excluded. 

This conclusion is corroborated by the results for the 
CRSOS model with AH max = 2, although the accuracy 
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of the data was poorer. In Fig. 7 we show zil^\ versus 
1/L for that model, with To also calculated using four 
different values of k in Eq. (10). 

Our results for the noise-reduced DT model do not 
provide useful information on dynamical exponents, sim- 
ilarly to the case of the roughness exponents. 



To. Consequently, we were not able to calculate accurate 
dynamical exponents in the two-dimensional case. 
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FIG. 7. Effective dynamical exponents Z(l,4,) versus 1/L for 



the 1 + 1- dimensional CRSOS model with AH„ 



2. Error 



bars (not shown) are smaller than Az 
for the largest L). 
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(a) 



(b) 



FIG. 8. Effective roughness exponents (a) o.(l,2) (obtained 
from the interface width) and (b) 2 ) (obtained from Wa) 
versus 1/L for the 2 + 1-dimensional CRSOS model with 
AH max = 1. Error bars are shown only when they are larger 
than the size of the data points. 



V. UNIVERSALITY OF AMPLITUDE RATIOS 



IV. SCALING EXPONENTS IN 
TWO-DIMENSIONAL SUBSTRATES 

In Figs. 8a and 8b we show Q(l,2) (Eq- 12) and a^ 2 ) 
(Eq. 16) for the two-dimensional CRSOS model with 
AH max = 1. Both linear fits give a = 0.662, which is 
very near the one-loop renormalization value a — 2/3 of 
the VLDS theory. Accounting for the error bars, which 
are particularly large for L = 256, we are not able to 
determine whether a — 2/3 is exact or not. On the 
other hand, confirming other authors' results [10], any 
difference from that value is probably smaller than the 
two-loops correction of Janssen [7], which is Aa « 0.014. 

Similarly to the one-dimensional case, the error bars 
of the data for the model with AH max — 2 are larger. 
Consequently, no discrepancy from the one-loop expo- 
nents could be detected too. 

The characteristic times To for the model with 
AH max = 1 were obtained in lattices with 16 < L < 128, 
but their values for the smallest lattices (L = 16 and 
L = 32) are very small, sometimes below tq = 1 (one 
monolayer). For L = 256, the accuracy of the interface 
widths data is not enough to provide reliable estimates of 



Evidence on the power-counting property of the mo- 
ments W n of the heights distribution of VLDS models 
was given in Sec. Ill by the estimates of a obtained from 
W2 and W4 . A clearer evidence is given here by the finite 
asymptotic estimates of the skewness and the kurtosis at 
the steady states. 

First we consider the models in 1 + 1 dimensions. 

In Figs. 9a and 9b we show the steady state skewness 
versus 1/L 1 / 2 for the CRSOS models with AH max = 
1 and AH max — 2, respectively. Except for the data 
for L = 1024, which have relatively large error bars, all 
points fall in almost perfect straight lines, which give the 
asymptotic value S = 0.32 ± 0.02 for both models. 

In Figs. 9c and 9d we show the steady state kurtosis 
versus 1/L 1 / 2 for the CRSOS models with AH max = 1 
and AH max = 2, respectively. Only the data for L < 
512 were shown because the error bars are much larger 
for L = 1024, not giving additional information on the 
evolution of Q. Reasonable linear fits are obtained with 
the last four data points in each case. The asymptotic 
estimate is Q = —0.11 ± 0.02 for both models. 

Our results for the competitive model (RD and CR- 
SOS) introduced in Sec. II also suggest that those ampli- 
tude ratios are universal for VLDS models. In that case, 
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there is no constraint on the difference of the heights of 
neighboring columns, but only a trend to suppress large 
heights differences. The coefficients ^4 and A4 in the 
corresponding continuous equation (Eq. 1) are probably 
different from those in the pure model (p = 0), as ob- 
tained in related competitive models [28,29]. In Figs. 
10a and 10b we show, respectively, S (L, t — > 00) and 
Q(L,t — > 00) as a function of 1/L 1 / 2 for the competitive 
model. The asymptotic estimates are S = 0.32 ± 0.02 
and Q s=s —0.1, which are near the previous estimates for 
the pure CRSOS model. 
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FIG. 9. Steady state skewness for the 1 + 1-dimensional 
CRSOS model with (a) AH max = 1 and (b)AH max = 2, and 
steady state kurtosis for that model with (c) AH max = 1 and 
(d) AHmax — 2, as functions of 1/L 1 / 2 . Dotted lines are least 
squares fits of the data. Error bars are shown only when they 
are larger than the size of the data points. 



In Figs. 10c and lOd we show, respectively, 
S (L,t — > 00) and Q{L,t—* 00) as a function of 1/L 1 / 2 
for the noise-reduced DT model in d = 1. There are 
several reasons for the large error bars of the kurtosis, 
particularly in the largest lattices. Firstly, as justified in 
Sec. Ill, fluctuations in the data for the DT model are 
typically large. Secondly, the relative fluctuations of the 
moments W n (Eq. 6) rapidly increase with the order n. 
Finally, while the size of the error bar of the kurtosis is 
the same of W4/W2 2 , the relative error significantly in- 
creases when the constant 3 is subtracted (Eq. 9). The 
relatively large errors in Figs. 9c and 9d (CRSOS mod- 
els) can also be explained along these lines. 
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FIG. 10. (a), (b): steady state skewness and kurtosis, re- 
spectively, as a function of 1/L 1 '' 2 , for the competitive model 
(CRSOS with AH max = 1 and RD); (c), (d): steady state 
skewness and kurtosis, respectively, as a function of 1/L 1 / 2 , 
for the DT model. Dotted lines are least squares fits of the 
data. Error bars are shown only when they are larger than 
the size of the data points. 



The trends of the data for the DT model in Figs. 10c 
and lOd are completely different from those of the CR- 
SOS models. We cannot exclude the possibility that the 
universality of the amplitude ratios be a special feature 
of CRSOS models and some simple extensions, like the 
above competitive model. However, the behavior of the 
scaling exponents of the DT model is also unusual, with 
no possible extrapolation to the expected region of the 
VLDS theory (a < 1, z < 3), as discussed in Sec. III. 
Consequently, the present results for the DT model, al- 
though not confirming the universality of the amplitude 
ratios, are not reliable to discard that hypothesis (the 
negative sign of the skewness is not a problem, since its 
sign changes with A4 - see related discussion in Ref. [13]). 

Now we turn to the CRSOS models in 2+1 dimensions. 

In Figs. 11a and lib we show the steady state skewness 
versus 1/L 1 / 2 for the CRSOS models with AH max = 1 
and AH max — 2, respectively. The asymptotic estimates 
are S = 0.19 ± 0.02 and S = 0.20 ± 0.02, which also sug- 
gest the universality of this quantity. In Figs. 11c and 
lid we show the steady state kurtosis versus 1/L 1 / 2 for 



the CRSOS models with AH r 



1 and AH r 



respectively. The asymptotic value Q — 0, which is the 
Gaussian value, is consistent with the error bars. Thus, 
in 2 + 1 dimensions, we also obtain evidence of univer- 
sality of the amplitude ratios for CRSOS models, which 
suggests this possibility for the whole VLDS class. 
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vergence of the DT data to the VLDS behavior. The hy- 
pothesis of a slow crossover is supported by the fact that 
the estimates of a for the DT model are significantly 
larger than the values predicted theoretically and con- 
firmed numerically (a < 1 in d = 1). Another possibility 
is that both CRSOS models and the competitive model 
have continuum representations with suitable combina- 
tions of coefficients which lead to the same forms of the 
heights distributions. 

We believe that the results of this work will motivate 
further studies, numerical and analytical, of the VLDS 
equation and related discrete models. The estimates of 
scaling exponents in d = 1 and the apparent universal- 
ity of amplitude ratios are some of the results that may 
eventually help one to validate approximations in ana- 
lytical works. On the other hand, numerical solutions of 
the VLDS equation or simulations of new discrete mod- 
els in this class would be relevant to broaden the present 
discussion. 



FIG. 11. Steady state skewness for the 2 + 1-dimensional 
CRSOS model with (a) AH max = 1 and (b)AH max = 2, and 
steady state kurtosis for that model with (c) AH ma x = 1 and 
(d) AH m ax = 2, as functions of 1/L 1 '' 2 . Dotted lines are least 
squares fits of the data. Error bars are shown only when they 
are larger than the size of the data points. 
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VI. SUMMARY AND CONCLUSION 



We studied numerically discrete growth models which 
belong to the VLDS class in 1 + 1 and 2 + 1 dimen- 
sions. Scaling exponents and steady state values of the 
skewness and the kurtosis, which characterize the heights 
distribution, were determined for those models. 

Results for the CRSOS model with AH max = 1 gave 
the roughness exponent a — 0.94 ± 0.02 and the dynam- 
ical exponent z = 2.88 ± 0.04 in d = 1. These estimates 
confirm the proposal of Janssen [7] that the exponents 
of the VLDS theory obtained from one-loop renormaliza- 
tion (a = 1 and z = 3) arc not exact. The corrections 
from two-loops calculations give a s=s 0.97 and z = 2.94, 
but they are obtained from expansions in 4 — d, which 
are not expected to provide accurate results for small d. 
On the other hand, the negative sign of the correction to 
one-loop results is consistent with our findings. In d = 2, 
our results are not able to exclude the one- loop values, 
confirming other authors' conclusions [10]. 

The estimates of the steady state skewness and kurtosis 
of the CRSOS models with AH max = 1 and AH max = 2 
and of the competitive model (RD versus CRSOS with 
AH max = 1) suggest that those amplitude ratios are uni- 
versal in the VLDS class. However, for the DT model in 
d = 1, which belongs to the same class, those quantities 
are very different from the suggested universal values. 
One possible reason for this discrepancy is the slow con- 
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